This script is used for summarizing IGC Simulation result (with tract length)
Read in data first
rm(list=ls()) # clean up workspace
IGC.path <- "/Users/xji3/FromCluster_IGCCodonSimulation01262016/"
paml.path <- "/Users/xji3/GitFolders/IGCCodonSimulation/"
IGC.geo.list <- c(3.0, 10.0, 50.0, 100.0, 500.0)
#IGC.geo.list <- c(1.0)
num.sim = 100
paralog = "YDR418W_YEL054C"
# Now read in PAML analysis of these simulated data set
data.path <- paste(paralog, "",sep = "/")
for (IGC.geo in IGC.geo.list){
summary_mat <- NULL
IGC.geo.path <- paste("IGCgeo_", toString(IGC.geo), ".0/", sep = "")
file.name <- paste("geo", paste(toString(IGC.geo), ".0", sep = ""), "paml", "unrooted", "summary.txt", sep = "_")
for (sim.num in 0:(num.sim - 1)){
summary_file <- paste(paml.path, file.name, sep = "")
if (file.exists(summary_file)){
all <- readLines(summary_file, n = -1)
col.names <- strsplit(all[1], ' ')[[1]][-1]
row.names <- strsplit(all[length(all)], ' ')[[1]][-1]
summary_mat <- as.matrix(read.table(summary_file,
row.names = row.names,
col.names = col.names))
}
}
assign(paste("PAML", paste(toString(IGC.geo), ".0", sep = ""), "summary", sep = "_"), summary_mat)
}
#Read.summary <- function(IGC.path, IGC.x3.path, IGC.geo.list, num.sim = 100, paralog = "YDR418W_YEL054C"){
## Now read simulation result
sim.path <- paste("SimulationSummary", paralog, "",sep = "/")
for (IGC.geo in IGC.geo.list){
summary_mat <- NULL
IGC.geo.path <- paste("IGCgeo_", toString(IGC.geo), ".0/", sep = "")
file.prefix <- paste(paralog, "MG94_geo", paste(toString(IGC.geo), ".0", sep = ""), "Sim", sep = "_")
file.suffix <- "summary.txt"
for (sim.num in 0:(num.sim - 1)){
file.name <- paste(file.prefix, toString(sim.num), file.suffix, sep = "_")
summary_file <- paste(IGC.path, sim.path, IGC.geo.path, file.name, sep = "")
if (file.exists(summary_file)){
all <- readLines(summary_file, n = -1)
col.names <- paste("geo", paste(toString(IGC.geo), ".0", sep = ""), "sim", toString(sim.num), sep = "_")
row.names <- strsplit(all[length(all)], ' ')[[1]][-1]
summary_mat <- cbind(summary_mat, as.matrix(read.table(summary_file,
row.names = row.names,
col.names = col.names)))
}
}
assign(paste("geo", paste(toString(IGC.geo), ".0", sep = ""), "summary", sep = "_"), summary_mat)
}
# Now read in 'real' % changes due to IGC in the simulation
data.path <- paste(paralog, "",sep = "/")
for (IGC.geo in IGC.geo.list){
summary_mat <- NULL
IGC.geo.path <- paste("IGCgeo_", toString(IGC.geo), ".0/", sep = "")
file.prefix <- paste(paralog, "MG94_geo", paste(toString(IGC.geo), ".0", sep = ""), "Sim", sep = "_")
file.suffix <- "short.log"
for (sim.num in 0:(num.sim - 1)){
file.name <- paste(file.prefix, toString(sim.num), file.suffix, sep = "_")
file.name <- paste("sim", paste(toString(sim.num), "/", file.prefix, sep = ""),
toString(sim.num), file.suffix, sep = "_")
summary_file <- paste(IGC.path, data.path, IGC.geo.path, file.name, sep = "")
if (file.exists(summary_file)){
all <- readLines(summary_file, n = -1)
col.names <- paste("geo", paste(toString(IGC.geo), ".0", sep = ""), "sim", toString(sim.num), sep = "_")
row.names <- strsplit(all[length(all)], ' ')[[1]][-1]
summary_mat <- cbind(summary_mat, as.matrix(read.table(summary_file,
row.names = row.names,
col.names = col.names)))
}
}
assign(paste("pIGC", paste(toString(IGC.geo), ".0", sep = ""), "summary", sep = "_"), summary_mat)
}
Now analyze the simulation results
parameters used in generating datasets are:
pi_a = 0.2983654, pi_c = 0.2095729, pi_g = 0.1871536, pi_t = 0.3049081
kappa = 8.4043336
omega = 1.0
tau = 1.409408 (0.42)
Tau estimate summary
# geo_3.0
summary(geo_3.0_summary["tau", ])
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.815 1.200 1.290 1.360 1.550 1.950
sd(geo_3.0_summary["tau", ])
## [1] 0.2547
# geo_10.0
summary(geo_10.0_summary["tau", ])
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.539 1.200 1.450 1.470 1.670 2.650
sd(geo_10.0_summary["tau", ])
## [1] 0.4084
# geo_50.0
summary(geo_50.0_summary["tau", ])
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.261 1.020 1.380 1.490 1.930 3.890
sd(geo_50.0_summary["tau", ])
## [1] 0.67
# geo_100.0
summary(geo_100.0_summary["tau", ])
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.114 1.090 1.600 1.620 1.940 4.500
sd(geo_100.0_summary["tau", ])
## [1] 0.779
# geo_500.0
summary(geo_500.0_summary["tau", ])
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.00 0.78 1.78 1.98 2.91 7.14
sd(geo_500.0_summary["tau", ])
## [1] 1.52
sd.list <- c(sd(geo_3.0_summary["tau", ]),
sd(geo_10.0_summary["tau", ]),sd(geo_50.0_summary["tau", ]),
sd(geo_100.0_summary["tau", ]), sd(geo_500.0_summary["tau", ]))
plot(IGC.geo.list, sd.list)
abline(h = 0.42)
mean.list <- c(mean(geo_3.0_summary["tau", ]),
mean(geo_10.0_summary["tau", ]), mean(geo_50.0_summary["tau", ]),
mean(geo_100.0_summary["tau", ]), mean(geo_500.0_summary["tau", ]))
plot(IGC.geo.list, mean.list)
abline(h = 1.409408)
Now summarized % changes due to IGC in simulation estimates
It was estimated 0.2131 (0.02884)
percent.IGC.geo.3.0 <- colSums(geo_3.0_summary[34:45, ] + geo_3.0_summary[46:57, ]) / (colSums(geo_3.0_summary[34:45, ] + geo_3.0_summary[46:57, ] + geo_3.0_summary[58:69, ]))
summary(percent.IGC.geo.3.0)
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.151 0.194 0.209 0.210 0.227 0.265
sd(percent.IGC.geo.3.0)
## [1] 0.02406
percent.IGC.geo.10.0 <- colSums(geo_10.0_summary[34:45, ] + geo_10.0_summary[46:57, ]) / (colSums(geo_10.0_summary[34:45, ] + geo_10.0_summary[46:57, ] + geo_10.0_summary[58:69, ]))
summary(percent.IGC.geo.10.0)
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.133 0.193 0.215 0.215 0.243 0.287
sd(percent.IGC.geo.10.0)
## [1] 0.03418
percent.IGC.geo.50.0 <- colSums(geo_50.0_summary[34:45, ] + geo_50.0_summary[46:57, ]) / (colSums(geo_50.0_summary[34:45, ] + geo_50.0_summary[46:57, ] + geo_50.0_summary[58:69, ]))
summary(percent.IGC.geo.50.0)
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.0678 0.1780 0.2140 0.2130 0.2580 0.3750
sd(percent.IGC.geo.50.0)
## [1] 0.05983
percent.IGC.geo.100.0 <- colSums(geo_100.0_summary[34:45, ] + geo_100.0_summary[46:57, ]) / (colSums(geo_100.0_summary[34:45, ] + geo_100.0_summary[46:57, ] + geo_100.0_summary[58:69, ]))
summary(percent.IGC.geo.100.0)
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.0286 0.1690 0.2320 0.2230 0.2670 0.3880
sd(percent.IGC.geo.100.0)
## [1] 0.06909
percent.IGC.geo.500.0 <- colSums(geo_500.0_summary[34:45, ] + geo_500.0_summary[46:57, ]) / (colSums(geo_500.0_summary[34:45, ] + geo_500.0_summary[46:57, ] + geo_500.0_summary[58:69, ]))
summary(percent.IGC.geo.500.0)
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.000 0.137 0.232 0.226 0.326 0.488
sd(percent.IGC.geo.500.0)
## [1] 0.1272
# plot mean % changes due to IGC v.s. IGC_geo
mean.percent.IGC.list <- c( mean(percent.IGC.geo.3.0),
mean(percent.IGC.geo.10.0), mean(percent.IGC.geo.50.0),
mean(percent.IGC.geo.100.0), mean(percent.IGC.geo.500.0))
plot(IGC.geo.list, mean.percent.IGC.list)
abline(h = 0.2131)
sd.percent.IGC.list <- c(sd(percent.IGC.geo.3.0),
sd(percent.IGC.geo.10.0), sd(percent.IGC.geo.50.0),
sd(percent.IGC.geo.100.0), sd(percent.IGC.geo.500.0))
plot(IGC.geo.list, sd.percent.IGC.list)
abline(h = 0.02884)
Now get histograms of tau estimate
# geo_3.0
hist(geo_3.0_summary["tau",])
# geo_10.0
hist(geo_10.0_summary["tau",])
# geo_50.0
hist(geo_50.0_summary["tau",])
# geo_100.0
hist(geo_100.0_summary["tau",])
# geo_500.0
hist(geo_500.0_summary["tau",])
Now get histograms of % change from IGC from estimate
# geo_3.0
hist(percent.IGC.geo.3.0)
# geo_10.0
hist(percent.IGC.geo.10.0)
# geo_50.0
hist(percent.IGC.geo.50.0)
# geo_100.0
hist(percent.IGC.geo.100.0)
# geo_500.0
hist(percent.IGC.geo.500.0)
Now summarize real % change due to IGC in simulation
It was estimated 0.2131 (0.02884) in real data set
pIGC.real.geo.3.0 <- pIGC_3.0_summary[1, ]
summary(pIGC.real.geo.3.0)
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.168 0.200 0.213 0.215 0.228 0.262
sd(pIGC.real.geo.3.0)
## [1] 0.01977
pIGC.real.geo.10.0 <- pIGC_10.0_summary[1, ]
summary(pIGC.real.geo.10.0)
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.139 0.198 0.214 0.215 0.238 0.279
sd(pIGC.real.geo.10.0)
## [1] 0.02781
pIGC.real.geo.50.0 <- pIGC_50.0_summary[1, ]
summary(pIGC.real.geo.50.0)
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.0698 0.1790 0.2140 0.2120 0.2520 0.3150
sd(pIGC.real.geo.50.0)
## [1] 0.0522
pIGC.real.geo.100.0 <- pIGC_100.0_summary[1, ]
summary(pIGC.real.geo.100.0)
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.0337 0.1790 0.2170 0.2170 0.2530 0.3580
sd(pIGC.real.geo.100.0)
## [1] 0.06202
pIGC.real.geo.500.0 <- pIGC_500.0_summary[1, ]
summary(pIGC.real.geo.500.0)
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.000 0.146 0.214 0.210 0.273 0.436
sd(pIGC.real.geo.500.0)
## [1] 0.1049
# plot mean % changes due to IGC v.s. IGC_geo
mean.pIGC.real.list <- c( mean(pIGC.real.geo.3.0),
mean(pIGC.real.geo.10.0), mean(pIGC.real.geo.50.0),
mean(pIGC.real.geo.100.0), mean(pIGC.real.geo.500.0))
plot(IGC.geo.list, mean.pIGC.real.list)
abline(h = 0.2131)
sd.pIGC.real.list <- c(sd(pIGC.real.geo.3.0),
sd(pIGC.real.geo.10.0), sd(pIGC.real.geo.50.0),
sd(pIGC.real.geo.100.0), sd(pIGC.real.geo.500.0))
plot(IGC.geo.list, sd.pIGC.real.list)
abline(h = 0.02884)
Now get histograms of real % change due to IGC from simulation
# geo_3.0
hist(pIGC.real.geo.3.0)
# geo_10.0
hist(pIGC.real.geo.10.0)
# geo_50.0
hist(pIGC.real.geo.50.0)
# geo_100.0
hist(pIGC.real.geo.100.0)
# geo_500.0
hist(pIGC.real.geo.500.0)
Now summarize difference between real % change due to IGC with estimated %
# geo_3.0
diff.pIGC.geo.3.0 <- percent.IGC.geo.3.0 - pIGC.real.geo.3.0
summary(diff.pIGC.geo.3.0)
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## -0.03820 -0.01680 -0.00743 -0.00451 0.00837 0.03130
sd(diff.pIGC.geo.3.0)
## [1] 0.01688
# geo_10.0
diff.pIGC.geo.10.0 <- percent.IGC.geo.10.0 - pIGC.real.geo.10.0
summary(diff.pIGC.geo.10.0)
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## -0.05680 -0.01590 0.00092 -0.00011 0.01410 0.05400
sd(diff.pIGC.geo.10.0)
## [1] 0.02105
# geo_50.0
diff.pIGC.geo.50.0 <- percent.IGC.geo.50.0 - pIGC.real.geo.50.0
summary(diff.pIGC.geo.50.0)
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## -0.05880 -0.01320 -0.00050 0.00039 0.01450 0.07750
sd(diff.pIGC.geo.50.0)
## [1] 0.02502
# geo_100.0
diff.pIGC.geo.100.0 <- percent.IGC.geo.100.0 - pIGC.real.geo.100.0
summary(diff.pIGC.geo.100.0)
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## -0.04790 -0.01170 0.00460 0.00546 0.01700 0.08940
sd(diff.pIGC.geo.100.0)
## [1] 0.02377
# geo_500.0
diff.pIGC.geo.500.0 <- percent.IGC.geo.500.0 - pIGC.real.geo.500.0
summary(diff.pIGC.geo.500.0)
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## -0.06150 -0.01080 0.00762 0.01600 0.04040 0.11700
sd(diff.pIGC.geo.500.0)
## [1] 0.03958
# plot mean % changes due to IGC v.s. IGC_geo
mean.diff.pIGC.list <- c(mean(diff.pIGC.geo.3.0),
mean(diff.pIGC.geo.10.0), mean(diff.pIGC.geo.50.0),
mean(diff.pIGC.geo.100.0), mean(diff.pIGC.geo.500.0))
plot(IGC.geo.list, mean.diff.pIGC.list)
abline(h = 0.0)
sd.diff.pIGC.list <- c(sd(diff.pIGC.geo.3.0),
sd(diff.pIGC.geo.10.0), sd(diff.pIGC.geo.50.0),
sd(diff.pIGC.geo.100.0), sd(diff.pIGC.geo.500.0))
plot(IGC.geo.list, sd.diff.pIGC.list)
abline(h = 0.02884)
Now get histograms of difference of % IGC from ‘real’ and estimate
# geo_3.0
hist(diff.pIGC.geo.3.0)
# geo_10.0
hist(diff.pIGC.geo.10.0)
# geo_50.0
hist(diff.pIGC.geo.50.0)
# geo_100.0
hist(diff.pIGC.geo.100.0)
# geo_500.0
hist(diff.pIGC.geo.500.0)
Now added PAML estimate of all simulated datasets. Start to show more comparisons.
Omega estimate from PAML
# geo_3.0
summary(PAML_3.0_summary["omega", ])
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.677 0.885 1.020 1.010 1.110 1.290
sd(PAML_3.0_summary["omega", ])
## [1] 0.138
# geo_10.0
summary(PAML_10.0_summary["omega", ])
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.605 0.914 0.971 0.996 1.100 1.480
sd(PAML_10.0_summary["omega", ])
## [1] 0.1516
# geo_50.0
summary(PAML_50.0_summary["omega", ])
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.725 0.933 1.020 1.020 1.120 1.490
sd(PAML_50.0_summary["omega", ])
## [1] 0.1596
# geo_100.0
summary(PAML_100.0_summary["omega", ])
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.688 0.904 0.987 0.996 1.060 1.500
sd(PAML_100.0_summary["omega", ])
## [1] 0.1447
# geo_500.0
summary(PAML_500.0_summary["omega", ])
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.606 0.896 1.010 1.020 1.120 1.510
sd(PAML_500.0_summary["omega", ])
## [1] 0.1677
PAML.omega.mean <- c(mean(PAML_3.0_summary["omega", ]), mean(PAML_10.0_summary["omega", ]),
mean(PAML_50.0_summary["omega", ]), mean(PAML_100.0_summary["omega", ]),
mean(PAML_500.0_summary["omega", ]))
PAML.omega.sd <- c(sd(PAML_3.0_summary["omega", ]), sd(PAML_10.0_summary["omega", ]),
sd(PAML_50.0_summary["omega", ]), sd(PAML_100.0_summary["omega", ]),
sd(PAML_500.0_summary["omega", ]))
plot(IGC.geo.list, PAML.omega.mean)
abline(h = 1.0)
plot(IGC.geo.list, PAML.omega.sd)
Omega estimate from IGC + MG94
# geo_3.0
summary(geo_3.0_summary["omega", ])
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.701 0.921 1.000 1.010 1.090 1.340
sd(geo_3.0_summary["omega", ])
## [1] 0.1234
# geo_10.0
summary(geo_10.0_summary["omega", ])
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.730 0.917 0.984 1.000 1.080 1.430
sd(geo_10.0_summary["omega", ])
## [1] 0.1407
# geo_50.0
summary(geo_50.0_summary["omega", ])
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.748 0.916 1.010 1.020 1.110 1.380
sd(geo_50.0_summary["omega", ])
## [1] 0.1436
# geo_100.0
summary(geo_100.0_summary["omega", ])
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.754 0.914 1.010 1.000 1.060 1.420
sd(geo_100.0_summary["omega", ])
## [1] 0.1314
# geo_500.0
summary(geo_500.0_summary["omega", ])
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.675 0.922 1.010 1.020 1.110 1.460
sd(geo_500.0_summary["omega", ])
## [1] 0.1481
geo.omega.mean <- c(mean(geo_3.0_summary["omega", ]), mean(geo_10.0_summary["omega", ]),
mean(geo_50.0_summary["omega", ]), mean(geo_100.0_summary["omega", ]),
mean(geo_500.0_summary["omega", ]))
geo.omega.sd <- c(sd(geo_3.0_summary["omega", ]), sd(geo_10.0_summary["omega", ]),
sd(geo_50.0_summary["omega", ]), sd(geo_100.0_summary["omega", ]),
sd(geo_500.0_summary["omega", ]))
plot(IGC.geo.list, geo.omega.mean)
abline(h = 1.0)
plot(IGC.geo.list, geo.omega.sd)
kappa estimate from PAML
# geo_3.0
summary(PAML_3.0_summary["kappa", ])
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 5.63 7.42 8.44 8.49 9.44 12.80
sd(PAML_3.0_summary["kappa", ])
## [1] 1.373
# geo_10.0
summary(PAML_10.0_summary["kappa", ])
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 6.10 7.58 8.11 8.36 9.03 14.30
sd(PAML_10.0_summary["kappa", ])
## [1] 1.39
# geo_50.0
summary(PAML_50.0_summary["kappa", ])
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 5.85 7.41 8.26 8.55 9.28 13.00
sd(PAML_50.0_summary["kappa", ])
## [1] 1.556
# geo_100.0
summary(PAML_100.0_summary["kappa", ])
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 5.29 7.33 8.30 8.49 9.53 13.00
sd(PAML_100.0_summary["kappa", ])
## [1] 1.572
# geo_500.0
summary(PAML_500.0_summary["kappa", ])
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 5.50 7.56 8.44 8.67 9.53 16.00
sd(PAML_500.0_summary["kappa", ])
## [1] 1.703
PAML.kappa.mean <- c(mean(PAML_3.0_summary["kappa", ]), mean(PAML_10.0_summary["kappa", ]),
mean(PAML_50.0_summary["kappa", ]), mean(PAML_100.0_summary["kappa", ]),
mean(PAML_500.0_summary["kappa", ]))
PAML.kappa.sd <- c(sd(PAML_3.0_summary["kappa", ]), sd(PAML_10.0_summary["kappa", ]),
sd(PAML_50.0_summary["kappa", ]), sd(PAML_100.0_summary["kappa", ]),
sd(PAML_500.0_summary["kappa", ]))
plot(IGC.geo.list, PAML.kappa.mean)
abline(h = 8.4043336)
plot(IGC.geo.list, PAML.kappa.sd)
kappa estimate from IGC + MG94
# geo_3.0
summary(geo_3.0_summary["kappa", ])
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 6.01 7.54 8.18 8.50 9.20 13.50
sd(geo_3.0_summary["kappa", ])
## [1] 1.321
# geo_10.0
summary(geo_10.0_summary["kappa", ])
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 5.61 7.59 8.21 8.43 8.99 13.10
sd(geo_10.0_summary["kappa", ])
## [1] 1.308
# geo_50.0
summary(geo_50.0_summary["kappa", ])
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 6.25 7.70 8.43 8.65 9.46 13.40
sd(geo_50.0_summary["kappa", ])
## [1] 1.44
# geo_100.0
summary(geo_100.0_summary["kappa", ])
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 5.42 7.51 8.43 8.59 9.50 13.00
sd(geo_100.0_summary["kappa", ])
## [1] 1.556
# geo_500.0
summary(geo_500.0_summary["kappa", ])
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 5.71 7.65 8.41 8.65 9.47 14.80
sd(geo_500.0_summary["kappa", ])
## [1] 1.613
geo.kappa.mean <- c(mean(geo_3.0_summary["kappa", ]), mean(geo_10.0_summary["kappa", ]),
mean(geo_50.0_summary["kappa", ]), mean(geo_100.0_summary["kappa", ]),
mean(geo_500.0_summary["kappa", ]))
geo.kappa.sd <- c(sd(geo_3.0_summary["kappa", ]), sd(geo_10.0_summary["kappa", ]),
sd(geo_50.0_summary["kappa", ]), sd(geo_100.0_summary["kappa", ]),
sd(geo_500.0_summary["kappa", ]))
plot(IGC.geo.list, geo.kappa.mean)
abline(h = 8.4043336)
plot(IGC.geo.list, geo.kappa.sd)
Now start to look at branch length estimates
Branch lengths used in simulation:
| Branch | blen |
|---|---|
| N0_N1: | 0.0197240946542 |
| N0_kluyveri | 0.215682181791 |
| N1_N2 | 0.20925129872 |
| N1_castellii | 0.171684721483 |
| N2_N3 | 0.0257112589202 |
| N2_bayanus | 0.0266075664688 |
| N3_N4 | 0.0321083243449 |
| N3_kudriavzevii | 0.0853588718458 |
| N4_N5 | 0.024947887926 |
| N4_mikatae | 0.0566627496729 |
| N5_cerevisiae | 0.0581451177847 |
| N5_paradoxus | 0.0218788166581 |
# N0_N1
IGC.N0.N1.mean.list <- c(mean(geo_3.0_summary["(N0,N1)", ]), mean(geo_10.0_summary["(N0,N1)", ]),
mean(geo_50.0_summary["(N0,N1)", ]), mean(geo_100.0_summary["(N0,N1)", ]),
mean(geo_500.0_summary["(N0,N1)", ]))
IGC.N0.N1.sd.list <- c(sd(geo_3.0_summary["(N0,N1)", ]), sd(geo_10.0_summary["(N0,N1)", ]),
sd(geo_50.0_summary["(N0,N1)", ]), sd(geo_100.0_summary["(N0,N1)", ]),
sd(geo_500.0_summary["(N0,N1)", ]))
PAML.N0.N1.paralog1.mean.list <- c(mean(PAML_3.0_summary["N0_N1", ]),
mean(PAML_10.0_summary["N0_N1", ]),
mean(PAML_50.0_summary["N0_N1", ]),
mean(PAML_100.0_summary["N0_N1", ]),
mean(PAML_500.0_summary["N0_N1", ]))
PAML.N0.N1.paralog1.sd.list <- c(sd(PAML_3.0_summary["N0_N1", ]),
sd(PAML_10.0_summary["N0_N1", ]),
sd(PAML_50.0_summary["N0_N1", ]),
sd(PAML_100.0_summary["N0_N1", ]),
sd(PAML_500.0_summary["N0_N1", ]))
PAML.N0.N1.paralog2.mean.list <- c(mean(PAML_3.0_summary["N0_N6", ]),
mean(PAML_10.0_summary["N0_N6", ]),
mean(PAML_50.0_summary["N0_N6", ]),
mean(PAML_100.0_summary["N0_N6", ]),
mean(PAML_500.0_summary["N0_N6", ]))
PAML.N0.N1.paralog2.sd.list <- c(sd(PAML_3.0_summary["N0_N6", ]),
sd(PAML_10.0_summary["N0_N6", ]),
sd(PAML_50.0_summary["N0_N6", ]),
sd(PAML_100.0_summary["N0_N6", ]),
sd(PAML_500.0_summary["N0_N6", ]))
# N0_kluyveri
IGC.N0.kluyveri.mean.list <- c(mean(geo_3.0_summary["(N0,kluyveri)", ]), mean(geo_10.0_summary["(N0,kluyveri)", ]),
mean(geo_50.0_summary["(N0,kluyveri)", ]), mean(geo_100.0_summary["(N0,kluyveri)", ]),
mean(geo_500.0_summary["(N0,kluyveri)", ]))
IGC.N0.kluyveri.sd.list <- c(sd(geo_3.0_summary["(N0,kluyveri)", ]), sd(geo_10.0_summary["(N0,kluyveri)", ]),
sd(geo_50.0_summary["(N0,kluyveri)", ]), sd(geo_100.0_summary["(N0,kluyveri)", ]),
sd(geo_500.0_summary["(N0,kluyveri)", ]))
PAML.N0.kluyveri.paralog1.mean.list <- c(mean(PAML_3.0_summary["N0_kluyveriYDR418W", ]),
mean(PAML_10.0_summary["N0_kluyveriYDR418W", ]),
mean(PAML_50.0_summary["N0_kluyveriYDR418W", ]),
mean(PAML_100.0_summary["N0_kluyveriYDR418W", ]),
mean(PAML_500.0_summary["N0_kluyveriYDR418W", ]))
PAML.N0.kluyveri.paralog1.sd.list <- c(sd(PAML_3.0_summary["N0_kluyveriYDR418W", ]),
sd(PAML_10.0_summary["N0_kluyveriYDR418W", ]),
sd(PAML_50.0_summary["N0_kluyveriYDR418W", ]),
sd(PAML_100.0_summary["N0_kluyveriYDR418W", ]),
sd(PAML_500.0_summary["N0_kluyveriYDR418W", ] ))
# N1_N2
IGC.N1.N2.mean.list <- c(mean(geo_3.0_summary["(N1,N2)", ]), mean(geo_10.0_summary["(N1,N2)", ]),
mean(geo_50.0_summary["(N1,N2)", ]), mean(geo_100.0_summary["(N1,N2)", ]),
mean(geo_500.0_summary["(N1,N2)", ]))
IGC.N1.N2.sd.list <- c(sd(geo_3.0_summary["(N1,N2)", ]), sd(geo_10.0_summary["(N1,N2)", ]),
sd(geo_50.0_summary["(N1,N2)", ]), sd(geo_100.0_summary["(N1,N2)", ]),
sd(geo_500.0_summary["(N1,N2)", ]))
PAML.N1.N2.paralog1.mean.list <- c(mean(PAML_3.0_summary["N1_N2", ]),
mean(PAML_10.0_summary["N1_N2", ]),
mean(PAML_50.0_summary["N1_N2", ]),
mean(PAML_100.0_summary["N1_N2", ]),
mean(PAML_500.0_summary["N1_N2", ]))
PAML.N1.N2.paralog1.sd.list <- c(sd(PAML_3.0_summary["N1_N2", ]),
sd(PAML_10.0_summary["N1_N2", ]),
sd(PAML_50.0_summary["N1_N2", ]),
sd(PAML_100.0_summary["N1_N2", ]),
sd(PAML_500.0_summary["N1_N2", ]))
PAML.N1.N2.paralog2.mean.list <- c(mean(PAML_3.0_summary["N6_N7", ]),
mean(PAML_10.0_summary["N6_N7", ]),
mean(PAML_50.0_summary["N6_N7", ]),
mean(PAML_100.0_summary["N6_N7", ]),
mean(PAML_500.0_summary["N6_N7", ]))
PAML.N1.N2.paralog2.sd.list <- c(sd(PAML_3.0_summary["N6_N7", ]),
sd(PAML_10.0_summary["N6_N7", ]),
sd(PAML_50.0_summary["N6_N7", ]),
sd(PAML_100.0_summary["N6_N7", ]),
sd(PAML_500.0_summary["N6_N7", ]))
# N1_castellii
IGC.N1.castellii.mean.list <- c(mean(geo_3.0_summary["(N1,castellii)", ]), mean(geo_10.0_summary["(N1,castellii)", ]),
mean(geo_50.0_summary["(N1,castellii)", ]), mean(geo_100.0_summary["(N1,castellii)", ]),
mean(geo_500.0_summary["(N1,castellii)", ]))
IGC.N1.castellii.sd.list <- c(sd(geo_3.0_summary["(N1,castellii)", ]), sd(geo_10.0_summary["(N1,castellii)", ]),
sd(geo_50.0_summary["(N1,castellii)", ]), sd(geo_100.0_summary["(N1,castellii)", ]),
sd(geo_500.0_summary["(N1,castellii)", ]))
PAML.N1.castellii.paralog1.mean.list <- c(mean(PAML_3.0_summary["N1_castelliiYDR418W", ]),
mean(PAML_10.0_summary["N1_castelliiYDR418W", ]),
mean(PAML_50.0_summary["N1_castelliiYDR418W", ]),
mean(PAML_100.0_summary["N1_castelliiYDR418W", ]),
mean(PAML_500.0_summary["N1_castelliiYDR418W", ]))
PAML.N1.castellii.paralog1.sd.list <- c(sd(PAML_3.0_summary["N1_castelliiYDR418W", ]),
sd(PAML_10.0_summary["N1_castelliiYDR418W", ]),
sd(PAML_50.0_summary["N1_castelliiYDR418W", ]),
sd(PAML_100.0_summary["N1_castelliiYDR418W", ]),
sd(PAML_500.0_summary["N1_castelliiYDR418W", ]))
PAML.N1.castellii.paralog2.mean.list <- c(mean(PAML_3.0_summary["N6_castelliiYEL054C", ]),
mean(PAML_10.0_summary["N6_castelliiYEL054C", ]),
mean(PAML_50.0_summary["N6_castelliiYEL054C", ]),
mean(PAML_100.0_summary["N6_castelliiYEL054C", ]),
mean(PAML_500.0_summary["N6_castelliiYEL054C", ]))
PAML.N1.castellii.paralog2.sd.list <- c(sd(PAML_3.0_summary["N6_castelliiYEL054C", ]),
sd(PAML_10.0_summary["N6_castelliiYEL054C", ]),
sd(PAML_50.0_summary["N6_castelliiYEL054C", ]),
sd(PAML_100.0_summary["N6_castelliiYEL054C", ]),
sd(PAML_500.0_summary["N6_castelliiYEL054C", ]))
# N2_N3
IGC.N2.N3.mean.list <- c(mean(geo_3.0_summary["(N2,N3)", ]), mean(geo_10.0_summary["(N2,N3)", ]),
mean(geo_50.0_summary["(N2,N3)", ]), mean(geo_100.0_summary["(N2,N3)", ]),
mean(geo_500.0_summary["(N2,N3)", ]))
IGC.N2.N3.sd.list <- c(sd(geo_3.0_summary["(N2,N3)", ]), sd(geo_10.0_summary["(N2,N3)", ]),
sd(geo_50.0_summary["(N2,N3)", ]), sd(geo_100.0_summary["(N2,N3)", ]),
sd(geo_500.0_summary["(N2,N3)", ]))
PAML.N2.N3.paralog1.mean.list <- c(mean(PAML_3.0_summary["N2_N3", ]),
mean(PAML_10.0_summary["N2_N3", ]),
mean(PAML_50.0_summary["N2_N3", ]),
mean(PAML_100.0_summary["N2_N3", ]),
mean(PAML_500.0_summary["N2_N3", ]))
PAML.N2.N3.paralog1.sd.list <- c(sd(PAML_3.0_summary["N2_N3", ]),
sd(PAML_10.0_summary["N2_N3", ]),
sd(PAML_50.0_summary["N2_N3", ]),
sd(PAML_100.0_summary["N2_N3", ]),
sd(PAML_500.0_summary["N2_N3", ]))
PAML.N2.N3.paralog2.mean.list <- c(mean(PAML_3.0_summary["N7_N8", ]),
mean(PAML_10.0_summary["N7_N8", ]),
mean(PAML_50.0_summary["N7_N8", ]),
mean(PAML_100.0_summary["N7_N8", ]),
mean(PAML_500.0_summary["N7_N8", ]))
PAML.N2.N3.paralog2.sd.list <- c(sd(PAML_3.0_summary["N7_N8", ]),
sd(PAML_10.0_summary["N7_N8", ]),
sd(PAML_50.0_summary["N7_N8", ]),
sd(PAML_100.0_summary["N7_N8", ]),
sd(PAML_500.0_summary["N7_N8", ]))
# N2_bayanus
IGC.N2.bayanus.mean.list <- c(mean(geo_3.0_summary["(N2,bayanus)", ]), mean(geo_10.0_summary["(N2,bayanus)", ]),
mean(geo_50.0_summary["(N2,bayanus)", ]), mean(geo_100.0_summary["(N2,bayanus)", ]),
mean(geo_500.0_summary["(N2,bayanus)", ]))
IGC.N2.bayanus.sd.list <- c(sd(geo_3.0_summary["(N2,bayanus)", ]), sd(geo_10.0_summary["(N2,bayanus)", ]),
sd(geo_50.0_summary["(N2,bayanus)", ]), sd(geo_100.0_summary["(N2,bayanus)", ]),
sd(geo_500.0_summary["(N2,bayanus)", ]))
PAML.N2.bayanus.paralog1.mean.list <- c(mean(PAML_3.0_summary["N2_bayanusYDR418W", ]),
mean(PAML_10.0_summary["N2_bayanusYDR418W", ]),
mean(PAML_50.0_summary["N2_bayanusYDR418W", ]),
mean(PAML_100.0_summary["N2_bayanusYDR418W", ]),
mean(PAML_500.0_summary["N2_bayanusYDR418W", ]))
PAML.N2.bayanus.paralog1.sd.list <- c(sd(PAML_3.0_summary["N2_bayanusYDR418W", ]),
sd(PAML_10.0_summary["N2_bayanusYDR418W", ]),
sd(PAML_50.0_summary["N2_bayanusYDR418W", ]),
sd(PAML_100.0_summary["N2_bayanusYDR418W", ]),
sd(PAML_500.0_summary["N2_bayanusYDR418W", ]))
PAML.N2.bayanus.paralog2.mean.list <- c(mean(PAML_3.0_summary["N7_bayanusYEL054C", ]),
mean(PAML_10.0_summary["N7_bayanusYEL054C", ]),
mean(PAML_50.0_summary["N7_bayanusYEL054C", ]),
mean(PAML_100.0_summary["N7_bayanusYEL054C", ]),
mean(PAML_500.0_summary["N7_bayanusYEL054C", ]))
PAML.N2.bayanus.paralog2.sd.list <- c(sd(PAML_3.0_summary["N7_bayanusYEL054C", ]),
sd(PAML_10.0_summary["N7_bayanusYEL054C", ]),
sd(PAML_50.0_summary["N7_bayanusYEL054C", ]),
sd(PAML_100.0_summary["N7_bayanusYEL054C", ]),
sd(PAML_500.0_summary["N7_bayanusYEL054C", ]))
# N3_N4
IGC.N3.N4.mean.list <- c(mean(geo_3.0_summary["(N3,N4)", ]), mean(geo_10.0_summary["(N3,N4)", ]),
mean(geo_50.0_summary["(N3,N4)", ]), mean(geo_100.0_summary["(N3,N4)", ]),
mean(geo_500.0_summary["(N3,N4)", ]))
IGC.N3.N4.sd.list <- c(sd(geo_3.0_summary["(N3,N4)", ]), sd(geo_10.0_summary["(N3,N4)", ]),
sd(geo_50.0_summary["(N3,N4)", ]), sd(geo_100.0_summary["(N3,N4)", ]),
sd(geo_500.0_summary["(N3,N4)", ]))
PAML.N3.N4.paralog1.mean.list <- c(mean(PAML_3.0_summary["N3_N4", ]),
mean(PAML_10.0_summary["N3_N4", ]),
mean(PAML_50.0_summary["N3_N4", ]),
mean(PAML_100.0_summary["N3_N4", ]),
mean(PAML_500.0_summary["N3_N4", ]))
PAML.N3.N4.paralog1.sd.list <- c(sd(PAML_3.0_summary["N3_N4", ]),
sd(PAML_10.0_summary["N3_N4", ]),
sd(PAML_50.0_summary["N3_N4", ]),
sd(PAML_100.0_summary["N3_N4", ]),
sd(PAML_500.0_summary["N3_N4", ]))
PAML.N3.N4.paralog2.mean.list <- c(mean(PAML_3.0_summary["N8_N9", ]),
mean(PAML_10.0_summary["N8_N9", ]),
mean(PAML_50.0_summary["N8_N9", ]),
mean(PAML_100.0_summary["N8_N9", ]),
mean(PAML_500.0_summary["N8_N9", ]))
PAML.N3.N4.paralog2.sd.list <- c(sd(PAML_3.0_summary["N8_N9", ]),
sd(PAML_10.0_summary["N8_N9", ]),
sd(PAML_50.0_summary["N8_N9", ]),
sd(PAML_100.0_summary["N8_N9", ]),
sd(PAML_500.0_summary["N8_N9", ]))
# N3_kudriavzevii
IGC.N3.kudriavzevii.mean.list <- c(mean(geo_3.0_summary["(N3,kudriavzevii)", ]), mean(geo_10.0_summary["(N3,kudriavzevii)", ]),
mean(geo_50.0_summary["(N3,kudriavzevii)", ]), mean(geo_100.0_summary["(N3,kudriavzevii)", ]),
mean(geo_500.0_summary["(N3,kudriavzevii)", ]))
IGC.N3.kudriavzevii.sd.list <- c(sd(geo_3.0_summary["(N3,kudriavzevii)", ]), sd(geo_10.0_summary["(N3,kudriavzevii)", ]),
sd(geo_50.0_summary["(N3,kudriavzevii)", ]), sd(geo_100.0_summary["(N3,kudriavzevii)", ]),
sd(geo_500.0_summary["(N3,kudriavzevii)", ]))
PAML.N3.kudriavzevii.paralog1.mean.list <- c(mean(PAML_3.0_summary["N3_kudriavzeviiYDR418W", ]),
mean(PAML_10.0_summary["N3_kudriavzeviiYDR418W", ]),
mean(PAML_50.0_summary["N3_kudriavzeviiYDR418W", ]),
mean(PAML_100.0_summary["N3_kudriavzeviiYDR418W", ]),
mean(PAML_500.0_summary["N3_kudriavzeviiYDR418W", ]))
PAML.N3.kudriavzevii.paralog1.sd.list <- c(sd(PAML_3.0_summary["N3_kudriavzeviiYDR418W", ]),
sd(PAML_10.0_summary["N3_kudriavzeviiYDR418W", ]),
sd(PAML_50.0_summary["N3_kudriavzeviiYDR418W", ]),
sd(PAML_100.0_summary["N3_kudriavzeviiYDR418W", ]),
sd(PAML_500.0_summary["N3_kudriavzeviiYDR418W", ]))
PAML.N3.kudriavzevii.paralog2.mean.list <- c(mean(PAML_3.0_summary["N8_kudriavzeviiYEL054C", ]),
mean(PAML_10.0_summary["N8_kudriavzeviiYEL054C", ]),
mean(PAML_50.0_summary["N8_kudriavzeviiYEL054C", ]),
mean(PAML_100.0_summary["N8_kudriavzeviiYEL054C", ]),
mean(PAML_500.0_summary["N8_kudriavzeviiYEL054C", ]))
PAML.N3.kudriavzevii.paralog2.sd.list <- c(sd(PAML_3.0_summary["N8_kudriavzeviiYEL054C", ]),
sd(PAML_10.0_summary["N8_kudriavzeviiYEL054C", ]),
sd(PAML_50.0_summary["N8_kudriavzeviiYEL054C", ]),
sd(PAML_100.0_summary["N8_kudriavzeviiYEL054C", ]),
sd(PAML_500.0_summary["N8_kudriavzeviiYEL054C", ]))
# N4_N5
IGC.N4.N5.mean.list <- c(mean(geo_3.0_summary["(N4,N5)", ]), mean(geo_10.0_summary["(N4,N5)", ]),
mean(geo_50.0_summary["(N4,N5)", ]), mean(geo_100.0_summary["(N4,N5)", ]),
mean(geo_500.0_summary["(N4,N5)", ]))
IGC.N4.N5.sd.list <- c(sd(geo_3.0_summary["(N4,N5)", ]), sd(geo_10.0_summary["(N4,N5)", ]),
sd(geo_50.0_summary["(N4,N5)", ]), sd(geo_100.0_summary["(N4,N5)", ]),
sd(geo_500.0_summary["(N4,N5)", ]))
PAML.N4.N5.paralog1.mean.list <- c(mean(PAML_3.0_summary["N4_N5", ]),
mean(PAML_10.0_summary["N4_N5", ]),
mean(PAML_50.0_summary["N4_N5", ]),
mean(PAML_100.0_summary["N4_N5", ]),
mean(PAML_500.0_summary["N4_N5", ]))
PAML.N4.N5.paralog1.sd.list <- c(sd(PAML_3.0_summary["N4_N5", ]),
sd(PAML_10.0_summary["N4_N5", ]),
sd(PAML_50.0_summary["N4_N5", ]),
sd(PAML_100.0_summary["N4_N5", ]),
sd(PAML_500.0_summary["N4_N5", ]))
PAML.N4.N5.paralog2.mean.list <- c(mean(PAML_3.0_summary["N9_N10", ]),
mean(PAML_10.0_summary["N9_N10", ]),
mean(PAML_50.0_summary["N9_N10", ]),
mean(PAML_100.0_summary["N9_N10", ]),
mean(PAML_500.0_summary["N9_N10", ]))
PAML.N4.N5.paralog2.sd.list <- c(sd(PAML_3.0_summary["N9_N10", ]),
sd(PAML_10.0_summary["N9_N10", ]),
sd(PAML_50.0_summary["N9_N10", ]),
sd(PAML_100.0_summary["N9_N10", ]),
sd(PAML_500.0_summary["N9_N10", ]))
# N4_mikatae
IGC.N4.mikatae.mean.list <- c(mean(geo_3.0_summary["(N4,mikatae)", ]), mean(geo_10.0_summary["(N4,mikatae)", ]),
mean(geo_50.0_summary["(N4,mikatae)", ]), mean(geo_100.0_summary["(N4,mikatae)", ]),
mean(geo_500.0_summary["(N4,mikatae)", ]))
IGC.N4.mikatae.sd.list <- c(sd(geo_3.0_summary["(N4,mikatae)", ]), sd(geo_10.0_summary["(N4,mikatae)", ]),
sd(geo_50.0_summary["(N4,mikatae)", ]), sd(geo_100.0_summary["(N4,mikatae)", ]),
sd(geo_500.0_summary["(N4,mikatae)", ]))
PAML.N4.mikatae.paralog1.mean.list <- c(mean(PAML_3.0_summary["N4_mikataeYDR418W", ]),
mean(PAML_10.0_summary["N4_mikataeYDR418W", ]),
mean(PAML_50.0_summary["N4_mikataeYDR418W", ]),
mean(PAML_100.0_summary["N4_mikataeYDR418W", ]),
mean(PAML_500.0_summary["N4_mikataeYDR418W", ]))
PAML.N4.mikatae.paralog1.sd.list <- c(sd(PAML_3.0_summary["N4_mikataeYDR418W", ]),
sd(PAML_10.0_summary["N4_mikataeYDR418W", ]),
sd(PAML_50.0_summary["N4_mikataeYDR418W", ]),
sd(PAML_100.0_summary["N4_mikataeYDR418W", ]),
sd(PAML_500.0_summary["N4_mikataeYDR418W", ]))
PAML.N4.mikatae.paralog2.mean.list <- c(mean(PAML_3.0_summary["N9_mikataeYEL054C", ]),
mean(PAML_10.0_summary["N9_mikataeYEL054C", ]),
mean(PAML_50.0_summary["N9_mikataeYEL054C", ]),
mean(PAML_100.0_summary["N9_mikataeYEL054C", ]),
mean(PAML_500.0_summary["N9_mikataeYEL054C", ]))
PAML.N4.mikatae.paralog2.sd.list <- c(sd(PAML_3.0_summary["N9_mikataeYEL054C", ]),
sd(PAML_10.0_summary["N9_mikataeYEL054C", ]),
sd(PAML_50.0_summary["N9_mikataeYEL054C", ]),
sd(PAML_100.0_summary["N9_mikataeYEL054C", ]),
sd(PAML_500.0_summary["N9_mikataeYEL054C", ]))
# N5_cerevisiae
IGC.N5.cerevisiae.mean.list <- c(mean(geo_3.0_summary["(N5,cerevisiae)", ]), mean(geo_10.0_summary["(N5,cerevisiae)", ]),
mean(geo_50.0_summary["(N5,cerevisiae)", ]), mean(geo_100.0_summary["(N5,cerevisiae)", ]),
mean(geo_500.0_summary["(N5,cerevisiae)", ]))
IGC.N5.cerevisiae.sd.list <- c(sd(geo_3.0_summary["(N5,cerevisiae)", ]), sd(geo_10.0_summary["(N5,cerevisiae)", ]),
sd(geo_50.0_summary["(N5,cerevisiae)", ]), sd(geo_100.0_summary["(N5,cerevisiae)", ]),
sd(geo_500.0_summary["(N5,cerevisiae)", ]))
PAML.N5.cerevisiae.paralog1.mean.list <- c(mean(PAML_3.0_summary["N5_cerevisiaeYDR418W", ]),
mean(PAML_10.0_summary["N5_cerevisiaeYDR418W", ]),
mean(PAML_50.0_summary["N5_cerevisiaeYDR418W", ]),
mean(PAML_100.0_summary["N5_cerevisiaeYDR418W", ]),
mean(PAML_500.0_summary["N5_cerevisiaeYDR418W", ]))
PAML.N5.cerevisiae.paralog1.sd.list <- c(sd(PAML_3.0_summary["N5_cerevisiaeYDR418W", ]),
sd(PAML_10.0_summary["N5_cerevisiaeYDR418W", ]),
sd(PAML_50.0_summary["N5_cerevisiaeYDR418W", ]),
sd(PAML_100.0_summary["N5_cerevisiaeYDR418W", ]),
sd(PAML_500.0_summary["N5_cerevisiaeYDR418W", ]))
PAML.N5.cerevisiae.paralog2.mean.list <- c(mean(PAML_3.0_summary["N10_cerevisiaeYEL054C", ]),
mean(PAML_10.0_summary["N10_cerevisiaeYEL054C", ]),
mean(PAML_50.0_summary["N10_cerevisiaeYEL054C", ]),
mean(PAML_100.0_summary["N10_cerevisiaeYEL054C", ]),
mean(PAML_500.0_summary["N10_cerevisiaeYEL054C", ]))
PAML.N5.cerevisiae.paralog2.sd.list <- c(sd(PAML_3.0_summary["N10_cerevisiaeYEL054C", ]),
sd(PAML_10.0_summary["N10_cerevisiaeYEL054C", ]),
sd(PAML_50.0_summary["N10_cerevisiaeYEL054C", ]),
sd(PAML_100.0_summary["N10_cerevisiaeYEL054C", ]),
sd(PAML_500.0_summary["N10_cerevisiaeYEL054C", ]))
# N5_paradoxus
IGC.N5.paradoxus.mean.list <- c(mean(geo_3.0_summary["(N5,paradoxus)", ]), mean(geo_10.0_summary["(N5,paradoxus)", ]),
mean(geo_50.0_summary["(N5,paradoxus)", ]), mean(geo_100.0_summary["(N5,paradoxus)", ]),
mean(geo_500.0_summary["(N5,paradoxus)", ]))
IGC.N5.paradoxus.sd.list <- c(sd(geo_3.0_summary["(N5,paradoxus)", ]), sd(geo_10.0_summary["(N5,paradoxus)", ]),
sd(geo_50.0_summary["(N5,paradoxus)", ]), sd(geo_100.0_summary["(N5,paradoxus)", ]),
sd(geo_500.0_summary["(N5,paradoxus)", ]))
PAML.N5.paradoxus.paralog1.mean.list <- c(mean(PAML_3.0_summary["N5_paradoxusYDR418W", ]),
mean(PAML_10.0_summary["N5_paradoxusYDR418W", ]),
mean(PAML_50.0_summary["N5_paradoxusYDR418W", ]),
mean(PAML_100.0_summary["N5_paradoxusYDR418W", ]),
mean(PAML_500.0_summary["N5_paradoxusYDR418W", ]))
PAML.N5.paradoxus.paralog1.sd.list <- c(sd(PAML_3.0_summary["N5_paradoxusYDR418W", ]),
sd(PAML_10.0_summary["N5_paradoxusYDR418W", ]),
sd(PAML_50.0_summary["N5_paradoxusYDR418W", ]),
sd(PAML_100.0_summary["N5_paradoxusYDR418W", ]),
sd(PAML_500.0_summary["N5_paradoxusYDR418W", ]))
PAML.N5.paradoxus.paralog2.mean.list <- c(mean(PAML_3.0_summary["N10_paradoxusYEL054C", ]),
mean(PAML_10.0_summary["N10_paradoxusYEL054C", ]),
mean(PAML_50.0_summary["N10_paradoxusYEL054C", ]),
mean(PAML_100.0_summary["N10_paradoxusYEL054C", ]),
mean(PAML_500.0_summary["N10_paradoxusYEL054C", ]))
PAML.N5.paradoxus.paralog2.sd.list <- c(sd(PAML_3.0_summary["N10_paradoxusYEL054C", ]),
sd(PAML_10.0_summary["N10_paradoxusYEL054C", ]),
sd(PAML_50.0_summary["N10_paradoxusYEL054C", ]),
sd(PAML_100.0_summary["N10_paradoxusYEL054C", ]),
sd(PAML_500.0_summary["N10_paradoxusYEL054C", ]))
Now plot all these tree branch length estimates from IGC + MG94 model
# N0_N1
plot(IGC.geo.list, IGC.N0.N1.mean.list)
abline(h = 0.0197240946542)
plot(IGC.geo.list, IGC.N0.N1.sd.list)
# N0_kluyveri
plot(IGC.geo.list, IGC.N0.kluyveri.mean.list)
abline(h = 0.215682181791)
plot(IGC.geo.list, IGC.N0.kluyveri.sd.list)
# N1_N2
plot(IGC.geo.list, IGC.N1.N2.mean.list)
abline(h = 0.20925129872)
plot(IGC.geo.list, IGC.N1.N2.sd.list)
# N1_castellii
plot(IGC.geo.list, IGC.N1.castellii.mean.list)
abline(h = 0.171684721483)
plot(IGC.geo.list, IGC.N1.castellii.sd.list)
# N2_N3
plot(IGC.geo.list, IGC.N2.N3.mean.list)
abline(h = 0.0257112589202)
plot(IGC.geo.list, IGC.N2.N3.sd.list)
# N2_bayanus
plot(IGC.geo.list, IGC.N2.bayanus.mean.list)
abline(h = 0.0266075664688)
plot(IGC.geo.list, IGC.N2.bayanus.sd.list)
# N3_N4
plot(IGC.geo.list, IGC.N3.N4.mean.list)
abline(h = 0.0321083243449)
plot(IGC.geo.list, IGC.N3.N4.sd.list)
# N3_kudriavzevii
plot(IGC.geo.list, IGC.N3.kudriavzevii.mean.list)
abline(h = 0.0853588718458)
plot(IGC.geo.list, IGC.N3.kudriavzevii.sd.list)
# N4_N5
plot(IGC.geo.list, IGC.N4.N5.mean.list)
abline(h = 0.024947887926)
plot(IGC.geo.list, IGC.N4.N5.sd.list)
# N4_mikatae
plot(IGC.geo.list, IGC.N4.mikatae.mean.list)
abline(h = 0.0566627496729)
plot(IGC.geo.list, IGC.N4.mikatae.sd.list)
# N5_cerevisiae
plot(IGC.geo.list, IGC.N5.cerevisiae.mean.list)
abline(h = 0.0581451177847)
plot(IGC.geo.list, IGC.N5.cerevisiae.sd.list)
# N5_paradoxus
plot(IGC.geo.list, IGC.N5.paradoxus.mean.list)
abline(h = 0.0218788166581)
plot(IGC.geo.list, IGC.N5.paradoxus.sd.list)
Summarize posterior #IGC + #mutation per site of each branch
post.branch.geo.3.0 <- (geo_3.0_summary[34:45, ] + geo_3.0_summary[46:57, ] + geo_3.0_summary[58:69, ])/geo_3.0_summary[1, ]
row.names(post.branch.geo.3.0) <- row.names(geo_3.0_summary)[10:21]
post.branch.geo.10.0 <- (geo_10.0_summary[34:45, ] + geo_10.0_summary[46:57, ] + geo_10.0_summary[58:69, ])/geo_10.0_summary[1, ]
row.names(post.branch.geo.10.0) <- row.names(geo_10.0_summary)[10:21]
post.branch.geo.50.0 <- (geo_50.0_summary[34:45, ] + geo_50.0_summary[46:57, ] + geo_50.0_summary[58:69, ])/geo_50.0_summary[1, ]
row.names(post.branch.geo.50.0) <- row.names(geo_50.0_summary)[10:21]
post.branch.geo.100.0 <- (geo_100.0_summary[34:45, ] + geo_100.0_summary[46:57, ] + geo_100.0_summary[58:69, ])/geo_100.0_summary[1, ]
row.names(post.branch.geo.100.0) <- row.names(geo_100.0_summary)[10:21]
post.branch.geo.500.0 <- (geo_500.0_summary[34:45, ] + geo_500.0_summary[46:57, ] + geo_500.0_summary[58:69, ])/geo_500.0_summary[1, ]
row.names(post.branch.geo.500.0) <- row.names(geo_500.0_summary)[10:21]
# N0_N1
mean.post.N0.N1.list <- c(mean(post.branch.geo.3.0["(N0,N1)", ]), mean(post.branch.geo.10.0["(N0,N1)", ]),
mean(post.branch.geo.50.0["(N0,N1)", ]), mean(post.branch.geo.100.0["(N0,N1)", ]),
mean(post.branch.geo.500.0["(N0,N1)", ])) / 2.0
# N0_kluyveri
mean.post.N0.kluyveri.list <- c(mean(post.branch.geo.3.0["(N0,kluyveri)", ]), mean(post.branch.geo.10.0["(N0,kluyveri)", ]),
mean(post.branch.geo.50.0["(N0,kluyveri)", ]), mean(post.branch.geo.100.0["(N0,kluyveri)", ]),
mean(post.branch.geo.500.0["(N0,kluyveri)", ]))
# N1_N2
mean.post.N1.N2.list <- c(mean(post.branch.geo.3.0["(N1,N2)", ]), mean(post.branch.geo.10.0["(N1,N2)", ]),
mean(post.branch.geo.50.0["(N1,N2)", ]), mean(post.branch.geo.100.0["(N1,N2)", ]),
mean(post.branch.geo.500.0["(N1,N2)", ])) / 2.0
# N1_castellii
mean.post.N1.castellii.list <- c(mean(post.branch.geo.3.0["(N1,castellii)", ]), mean(post.branch.geo.10.0["(N1,castellii)", ]),
mean(post.branch.geo.50.0["(N1,castellii)", ]), mean(post.branch.geo.100.0["(N1,castellii)", ]),
mean(post.branch.geo.500.0["(N1,castellii)", ])) / 2.0
# N2_N3
mean.post.N2.N3.list <- c(mean(post.branch.geo.3.0["(N2,N3)", ]), mean(post.branch.geo.10.0["(N2,N3)", ]),
mean(post.branch.geo.50.0["(N2,N3)", ]), mean(post.branch.geo.100.0["(N2,N3)", ]),
mean(post.branch.geo.500.0["(N2,N3)", ])) / 2.0
# N2_bayanus
mean.post.N2.bayanus.list <- c(mean(post.branch.geo.3.0["(N2,bayanus)", ]), mean(post.branch.geo.10.0["(N2,bayanus)", ]),
mean(post.branch.geo.50.0["(N2,bayanus)", ]), mean(post.branch.geo.100.0["(N2,bayanus)", ]),
mean(post.branch.geo.500.0["(N2,bayanus)", ])) / 2.0
# N3_N4
mean.post.N3.N4.list <- c(mean(post.branch.geo.3.0["(N3,N4)", ]), mean(post.branch.geo.10.0["(N3,N4)", ]),
mean(post.branch.geo.50.0["(N3,N4)", ]), mean(post.branch.geo.100.0["(N3,N4)", ]),
mean(post.branch.geo.500.0["(N3,N4)", ])) / 2.0
# N3_kudriavzevii
mean.post.N3.kudriavzevii.list <- c(mean(post.branch.geo.3.0["(N3,kudriavzevii)", ]), mean(post.branch.geo.10.0["(N3,kudriavzevii)", ]),
mean(post.branch.geo.50.0["(N3,kudriavzevii)", ]), mean(post.branch.geo.100.0["(N3,kudriavzevii)", ]),
mean(post.branch.geo.500.0["(N3,kudriavzevii)", ])) / 2.0
# N4_N5
mean.post.N4.N5.list <- c(mean(post.branch.geo.3.0["(N4,N5)", ]), mean(post.branch.geo.10.0["(N4,N5)", ]),
mean(post.branch.geo.50.0["(N4,N5)", ]), mean(post.branch.geo.100.0["(N4,N5)", ]),
mean(post.branch.geo.500.0["(N4,N5)", ])) / 2.0
# N4_mikatae
mean.post.N4.mikatae.list <- c(mean(post.branch.geo.3.0["(N4,mikatae)", ]), mean(post.branch.geo.10.0["(N4,mikatae)", ]),
mean(post.branch.geo.50.0["(N4,mikatae)", ]), mean(post.branch.geo.100.0["(N4,mikatae)", ]),
mean(post.branch.geo.500.0["(N4,mikatae)", ])) / 2.0
# N5_cerevisiae
mean.post.N5.cerevisiae.list <- c(mean(post.branch.geo.3.0["(N5,cerevisiae)", ]), mean(post.branch.geo.10.0["(N5,cerevisiae)", ]),
mean(post.branch.geo.50.0["(N5,cerevisiae)", ]), mean(post.branch.geo.100.0["(N5,cerevisiae)", ]),
mean(post.branch.geo.500.0["(N5,cerevisiae)", ])) / 2.0
# N5_paradoxus
mean.post.N5.paradoxus.list <- c(mean(post.branch.geo.3.0["(N5,paradoxus)", ]), mean(post.branch.geo.10.0["(N5,paradoxus)", ]),
mean(post.branch.geo.50.0["(N5,paradoxus)", ]), mean(post.branch.geo.100.0["(N5,paradoxus)", ]),
mean(post.branch.geo.500.0["(N5,paradoxus)", ])) / 2.0
# N0_N1
sd.post.N0.N1.list <- c(sd(post.branch.geo.3.0["(N0,N1)", ]), sd(post.branch.geo.10.0["(N0,N1)", ]),
sd(post.branch.geo.50.0["(N0,N1)", ]), sd(post.branch.geo.100.0["(N0,N1)", ]),
sd(post.branch.geo.500.0["(N0,N1)", ])) / 2.0
# N0_kluyveri
sd.post.N0.kluyveri.list <- c(sd(post.branch.geo.3.0["(N0,kluyveri)", ]), sd(post.branch.geo.10.0["(N0,kluyveri)", ]),
sd(post.branch.geo.50.0["(N0,kluyveri)", ]), sd(post.branch.geo.100.0["(N0,kluyveri)", ]),
sd(post.branch.geo.500.0["(N0,kluyveri)", ]))
# N1_N2
sd.post.N1.N2.list <- c(sd(post.branch.geo.3.0["(N1,N2)", ]), sd(post.branch.geo.10.0["(N1,N2)", ]),
sd(post.branch.geo.50.0["(N1,N2)", ]), sd(post.branch.geo.100.0["(N1,N2)", ]),
sd(post.branch.geo.500.0["(N1,N2)", ])) / 2.0
# N1_castellii
sd.post.N1.castellii.list <- c(sd(post.branch.geo.3.0["(N1,castellii)", ]), sd(post.branch.geo.10.0["(N1,castellii)", ]),
sd(post.branch.geo.50.0["(N1,castellii)", ]), sd(post.branch.geo.100.0["(N1,castellii)", ]),
sd(post.branch.geo.500.0["(N1,castellii)", ])) / 2.0
# N2_N3
sd.post.N2.N3.list <- c(sd(post.branch.geo.3.0["(N2,N3)", ]), sd(post.branch.geo.10.0["(N2,N3)", ]),
sd(post.branch.geo.50.0["(N2,N3)", ]), sd(post.branch.geo.100.0["(N2,N3)", ]),
sd(post.branch.geo.500.0["(N2,N3)", ])) / 2.0
# N2_bayanus
sd.post.N2.bayanus.list <- c(sd(post.branch.geo.3.0["(N2,bayanus)", ]), sd(post.branch.geo.10.0["(N2,bayanus)", ]),
sd(post.branch.geo.50.0["(N2,bayanus)", ]), sd(post.branch.geo.100.0["(N2,bayanus)", ]),
sd(post.branch.geo.500.0["(N2,bayanus)", ])) / 2.0
# N3_N4
sd.post.N3.N4.list <- c(sd(post.branch.geo.3.0["(N3,N4)", ]), sd(post.branch.geo.10.0["(N3,N4)", ]),
sd(post.branch.geo.50.0["(N3,N4)", ]), sd(post.branch.geo.100.0["(N3,N4)", ]),
sd(post.branch.geo.500.0["(N3,N4)", ])) / 2.0
# N3_kudriavzevii
sd.post.N3.kudriavzevii.list <- c(sd(post.branch.geo.3.0["(N3,kudriavzevii)", ]), sd(post.branch.geo.10.0["(N3,kudriavzevii)", ]),
sd(post.branch.geo.50.0["(N3,kudriavzevii)", ]), sd(post.branch.geo.100.0["(N3,kudriavzevii)", ]),
sd(post.branch.geo.500.0["(N3,kudriavzevii)", ])) / 2.0
# N4_N5
sd.post.N4.N5.list <- c(sd(post.branch.geo.3.0["(N4,N5)", ]), sd(post.branch.geo.10.0["(N4,N5)", ]),
sd(post.branch.geo.50.0["(N4,N5)", ]), sd(post.branch.geo.100.0["(N4,N5)", ]),
sd(post.branch.geo.500.0["(N4,N5)", ])) / 2.0
# N4_mikatae
sd.post.N4.mikatae.list <- c(sd(post.branch.geo.3.0["(N4,mikatae)", ]), sd(post.branch.geo.10.0["(N4,mikatae)", ]),
sd(post.branch.geo.50.0["(N4,mikatae)", ]), sd(post.branch.geo.100.0["(N4,mikatae)", ]),
sd(post.branch.geo.500.0["(N4,mikatae)", ])) / 2.0
# N5_cerevisiae
sd.post.N5.cerevisiae.list <- c(sd(post.branch.geo.3.0["(N5,cerevisiae)", ]), sd(post.branch.geo.10.0["(N5,cerevisiae)", ]),
sd(post.branch.geo.50.0["(N5,cerevisiae)", ]), sd(post.branch.geo.100.0["(N5,cerevisiae)", ]),
sd(post.branch.geo.500.0["(N5,cerevisiae)", ])) / 2.0
# N5_paradoxus
sd.post.N5.paradoxus.list <- c(sd(post.branch.geo.3.0["(N5,paradoxus)", ]), sd(post.branch.geo.10.0["(N5,paradoxus)", ]),
sd(post.branch.geo.50.0["(N5,paradoxus)", ]), sd(post.branch.geo.100.0["(N5,paradoxus)", ]),
sd(post.branch.geo.500.0["(N5,paradoxus)", ])) / 2.0
Now plot same branch length estimates from each paralog in paml results and posterior branch length from IGC + MG94 model
# N0_N1
matplot(IGC.geo.list, cbind(PAML.N0.N1.paralog1.mean.list, PAML.N0.N1.paralog2.mean.list, mean.post.N0.N1.list),
type = c("p"), pch = 1, col = 1:3, ylab = "mean N0.N1 estimate")
abline(h = 0.0197240946542)
legend("right", legend = c("paralog1", "paralog2", "posterior IGC"), col = 1:3, pch = 1)
matplot(IGC.geo.list, cbind(PAML.N0.N1.paralog1.sd.list, PAML.N0.N1.paralog2.sd.list, sd.post.N0.N1.list),
type = c("p"), pch = 1, col = 1:3, ylab = "sd N0.N1 estimate")
legend("right", legend = c("paralog1", "paralog2", "posterior IGC"), col = 1:3, pch = 1)
# N0_kluyveri
matplot(IGC.geo.list, cbind(PAML.N0.kluyveri.paralog1.mean.list, mean.post.N0.kluyveri.list),
type = c("p"), pch = 1, col = 1:3, ylab = "mean N0.kluyveri estimate")
abline(h = 0.215682181791)
legend("right", legend = c("paralog1", "posterior IGC"), col = 1:3, pch = 1)
matplot(IGC.geo.list, cbind(PAML.N0.kluyveri.paralog1.sd.list, sd.post.N0.kluyveri.list),
type = c("p"), pch = 1, col = 1:3, ylab = "sd N0.kluyveri estimate")
legend("right", legend = c("paralog1", "posterior IGC"), col = 1:3, pch = 1)
# N1_N2
matplot(IGC.geo.list, cbind(PAML.N1.N2.paralog1.mean.list, PAML.N1.N2.paralog2.mean.list, mean.post.N1.N2.list),
type = c("p"), pch = 1, col = 1:3, ylab = "mean N1.N2 estimate")
abline(h = 0.20925129872)
legend("right", legend = c("paralog1", "paralog2", "posterior IGC"), col = 1:3, pch = 1)
matplot(IGC.geo.list, cbind(PAML.N1.N2.paralog1.sd.list, PAML.N1.N2.paralog2.sd.list, sd.post.N1.N2.list),
type = c("p"), pch = 1, col = 1:3, ylab = "sd N1.N2 estimate")
legend("right", legend = c("paralog1", "paralog2", "posterior IGC"), col = 1:3, pch = 1)
# N1_castellii
matplot(IGC.geo.list, cbind(PAML.N1.castellii.paralog1.mean.list, PAML.N1.castellii.paralog2.mean.list, mean.post.N1.castellii.list),
type = c("p"), pch = 1, col = 1:3, ylab = "mean N1.castellii estimate")
abline(h = 0.171684721483)
legend("right", legend = c("paralog1", "paralog2", "posterior IGC"), col = 1:3, pch = 1)
matplot(IGC.geo.list, cbind(PAML.N1.castellii.paralog1.sd.list, PAML.N1.castellii.paralog2.sd.list, sd.post.N1.castellii.list),
type = c("p"), pch = 1, col = 1:3, ylab = "sd N1.castellii estimate")
legend("right", legend = c("paralog1", "paralog2", "posterior IGC"), col = 1:3, pch = 1)
# N2_N3
matplot(IGC.geo.list, cbind(PAML.N2.N3.paralog1.mean.list, PAML.N2.N3.paralog2.mean.list, mean.post.N2.N3.list),
type = c("p"), pch = 1, col = 1:3, ylab = "mean N2.N3 estimate")
abline(h = 0.0257112589202)
legend("right", legend = c("paralog1", "paralog2", "posterior IGC"), col = 1:3, pch = 1)
matplot(IGC.geo.list, cbind(PAML.N2.N3.paralog1.sd.list, PAML.N2.N3.paralog2.sd.list, sd.post.N2.N3.list),
type = c("p"), pch = 1, col = 1:3, ylab = "sd N2.N3 estimate")
legend("right", legend = c("paralog1", "paralog2", "posterior IGC"), col = 1:3, pch = 1)
# N2_bayanus
matplot(IGC.geo.list, cbind(PAML.N2.bayanus.paralog1.mean.list, PAML.N2.bayanus.paralog2.mean.list, mean.post.N2.bayanus.list),
type = c("p"), pch = 1, col = 1:3, ylab = "mean N2.bayanus estimate")
abline(h = 0.0266075664688)
legend("right", legend = c("paralog1", "paralog2", "posterior IGC"), col = 1:3, pch = 1)
matplot(IGC.geo.list, cbind(PAML.N2.bayanus.paralog1.sd.list, PAML.N2.bayanus.paralog2.sd.list, sd.post.N2.bayanus.list),
type = c("p"), pch = 1, col = 1:3, ylab = "sd N2.bayanus estimate")
legend("right", legend = c("paralog1", "paralog2", "posterior IGC"), col = 1:3, pch = 1)
# N3_N4
matplot(IGC.geo.list, cbind(PAML.N3.N4.paralog1.mean.list, PAML.N3.N4.paralog2.mean.list, mean.post.N3.N4.list),
type = c("p"), pch = 1, col = 1:3, ylab = "mean N3.N4 estimate")
abline(h = 0.0321083243449)
legend("right", legend = c("paralog1", "paralog2", "posterior IGC"), col = 1:3, pch = 1)
matplot(IGC.geo.list, cbind(PAML.N3.N4.paralog1.sd.list, PAML.N3.N4.paralog2.sd.list, sd.post.N3.N4.list),
type = c("p"), pch = 1, col = 1:3, ylab = "sd N3.N4 estimate")
legend("right", legend = c("paralog1", "paralog2", "posterior IGC"), col = 1:3, pch = 1)
# N3_kudriavzevii
matplot(IGC.geo.list, cbind(PAML.N3.kudriavzevii.paralog1.mean.list, PAML.N3.kudriavzevii.paralog2.mean.list, mean.post.N3.kudriavzevii.list),
type = c("p"), pch = 1, col = 1:3, ylab = "mean N3.kudriavzevii estimate")
abline(h = 0.0853588718458)
legend("right", legend = c("paralog1", "paralog2", "posterior IGC"), col = 1:3, pch = 1)
matplot(IGC.geo.list, cbind(PAML.N3.kudriavzevii.paralog1.sd.list, PAML.N3.kudriavzevii.paralog2.sd.list, sd.post.N3.kudriavzevii.list),
type = c("p"), pch = 1, col = 1:3, ylab = "sd N3.kudriavzevii estimate")
legend("right", legend = c("paralog1", "paralog2", "posterior IGC"), col = 1:3, pch = 1)
# N4_N5
matplot(IGC.geo.list, cbind(PAML.N4.N5.paralog1.mean.list, PAML.N4.N5.paralog2.mean.list, mean.post.N4.N5.list),
type = c("p"), pch = 1, col = 1:3, ylab = "mean N4.N5 estimate")
abline(h = 0.024947887926)
legend("right", legend = c("paralog1", "paralog2", "posterior IGC"), col = 1:3, pch = 1)
matplot(IGC.geo.list, cbind(PAML.N4.N5.paralog1.sd.list, PAML.N4.N5.paralog2.sd.list, sd.post.N4.N5.list),
type = c("p"), pch = 1, col = 1:3, ylab = "sd N4.N5 estimate")
legend("right", legend = c("paralog1", "paralog2", "posterior IGC"), col = 1:3, pch = 1)
# N4_mikatae
matplot(IGC.geo.list, cbind(PAML.N4.mikatae.paralog1.mean.list, PAML.N4.mikatae.paralog2.mean.list, mean.post.N4.mikatae.list),
type = c("p"), pch = 1, col = 1:3, ylab = "mean N4.mikatae estimate")
abline(h = 0.0566627496729)
legend("right", legend = c("paralog1", "paralog2", "posterior IGC"), col = 1:3, pch = 1)
matplot(IGC.geo.list, cbind(PAML.N4.mikatae.paralog1.sd.list, PAML.N4.mikatae.paralog2.sd.list, sd.post.N4.mikatae.list),
type = c("p"), pch = 1, col = 1:3, ylab = "sd N4.mikatae estimate")
legend("right", legend = c("paralog1", "paralog2", "posterior IGC"), col = 1:3, pch = 1)
# N5_cerevisiae
matplot(IGC.geo.list, cbind(PAML.N5.cerevisiae.paralog1.mean.list, PAML.N5.cerevisiae.paralog2.mean.list, mean.post.N5.cerevisiae.list),
type = c("p"), pch = 1, col = 1:3, ylab = "mean N5.cerevisiae estimate")
abline(h = 0.0581451177847)
legend("right", legend = c("paralog1", "paralog2", "posterior IGC"), col = 1:3, pch = 1)
matplot(IGC.geo.list, cbind(PAML.N5.cerevisiae.paralog1.sd.list, PAML.N5.cerevisiae.paralog2.sd.list, sd.post.N5.cerevisiae.list),
type = c("p"), pch = 1, col = 1:3, ylab = "sd N5.cerevisiae estimate")
legend("right", legend = c("paralog1", "paralog2", "posterior IGC"), col = 1:3, pch = 1)
# N5_paradoxus
matplot(IGC.geo.list, cbind(PAML.N5.paradoxus.paralog1.mean.list, PAML.N5.paradoxus.paralog2.mean.list, mean.post.N5.paradoxus.list),
type = c("p"), pch = 1, col = 1:3, ylab = "mean N5.paradoxus estimate")
abline(h = 0.0218788166581)
legend("right", legend = c("paralog1", "paralog2", "posterior IGC"), col = 1:3, pch = 1)
matplot(IGC.geo.list, cbind(PAML.N5.paradoxus.paralog1.sd.list, PAML.N5.paradoxus.paralog2.sd.list, sd.post.N5.paradoxus.list),
type = c("p"), pch = 1, col = 1:3, ylab = "sd N5.paradoxus estimate")
legend("right", legend = c("paralog1", "paralog2", "posterior IGC"), col = 1:3, pch = 1)
Now save workspace.
#save.image("/Users/xji3/GitFolders/IGCCodonSimulation/All.RData")